{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "70185b09",
   "metadata": {},
   "source": [
    "### Test code for SV window enrichment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3b713043",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np\n",
    "import pysam \n",
    "from pybedtools import BedTool\n",
    "from pathlib import Path\n",
    "from io import StringIO\n",
    "from functools import reduce"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "83256a69",
   "metadata": {},
   "outputs": [],
   "source": [
    "wkdir = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3\"\n",
    "wkdir_path = Path(wkdir)\n",
    "\n",
    "chrom=\"chr21\"\n",
    "vcf_path= \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/filtered_merged_gs_svp_10728.vcf.gz\"\n",
    "wgs_rna_paired_smpls_path= wkdir_path.joinpath(\"1_rna_seq_qc/wgs_rna_match/paired_wgs_rna_postqc_prioritise_wgs.tsv\")\n",
    "ge_matrix_flat_path = wkdir_path.joinpath(\"2_misexp_qc/misexp_gene_cov_corr/tpm_zscore_4568smpls_8610genes_flat_misexp_corr_qc.csv\")\n",
    "gencode_path = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/reference/gencode/gencode.v31.annotation.sorted.gtf.gz\"\n",
    "sv_info_path = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/info_table/final_sites_critical_info_allele.txt\"\n",
    "vep_msc_path = wkdir_path.joinpath(\"4_vrnt_enrich/sv_vep/msc/SV_vep_hg38_msc_parsed.tsv\")\n",
    "root_dir = wkdir_path.joinpath(\"4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "30c843b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "af_cutoff_list = [0, 0.01]\n",
    "af_lower, af_upper = af_cutoff_list\n",
    "window_start = 0\n",
    "window_step_size = 200000 \n",
    "window_max = 1000000 \n",
    "z_cutoff_list = [2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6b012d9e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inputs:\n",
      "- Chromosome: chr21\n",
      "- AF range: 0-0.01\n",
      "- Window start: 0\n",
      "- Window size: 200000\n",
      "- Max window: 1000000\n",
      "- Z-score cutoffs: [2]\n",
      "\n",
      "Gene IDs passing filters: 8610\n",
      "RNA-seq sample IDs passing QC: 4568\n",
      "\n",
      "Loading input VCF ...\n",
      "VCF loaded.\n",
      "Subset VCF to samples with RNA-seq ...\n",
      "Number of samples in VCF with RNA ID and passing QC: 2640\n",
      "Number of RNA IDs passing QC: 2640\n",
      "Number of genes pass QC on chr21: 148\n",
      "Sample IDs in gene expression matrix with SV calls: 2640\n",
      "SV types in SV info file: ['DEL', 'MEI', 'DUP', 'INV']\n"
     ]
    }
   ],
   "source": [
    "print(\"Inputs:\")\n",
    "print(f\"- Chromosome: {chrom}\")\n",
    "print(f\"- AF range: {af_lower}-{af_upper}\")\n",
    "print(f\"- Window start: {window_start}\")\n",
    "print(f\"- Window size: {window_step_size}\")\n",
    "print(f\"- Max window: {window_max}\")\n",
    "print(f\"- Z-score cutoffs: {z_cutoff_list}\")\n",
    "print(\"\")\n",
    "\n",
    "# constants \n",
    "tpm_cutoff = 0.5 \n",
    "\n",
    "# create root directory \n",
    "root_dir_path = Path(root_dir)\n",
    "root_dir_path.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "### Collect samples with VCF calls and RNA sequencing \n",
    "# read in gene expression file \n",
    "ge_matrix_flat_df = pd.read_csv(ge_matrix_flat_path, sep=\",\")\n",
    "gene_id_pass_qc_set = set(ge_matrix_flat_df.gene_id.unique())\n",
    "print(f\"Gene IDs passing filters: {len(gene_id_pass_qc_set)}\")\n",
    "smpl_id_pass_qc_set = set(ge_matrix_flat_df.rna_id.unique())\n",
    "print(f\"RNA-seq sample IDs passing QC: {len(smpl_id_pass_qc_set)}\")\n",
    "print(\"\")\n",
    "# egan ID, RNA ID sample links \n",
    "wgs_rna_paired_smpls_df = pd.read_csv(wgs_rna_paired_smpls_path, sep=\"\\t\")\n",
    "egan_ids_with_rna = wgs_rna_paired_smpls_df[wgs_rna_paired_smpls_df.rna_id.isin(smpl_id_pass_qc_set)].egan_id.tolist()\n",
    "# load VCF and subset to EGAN IDs with RNA \n",
    "print(\"Loading input VCF ...\")\n",
    "vcf_path = vcf_path\n",
    "vcf = pysam.VariantFile(vcf_path, mode = \"r\")\n",
    "print(\"VCF loaded.\")\n",
    "print(\"Subset VCF to samples with RNA-seq ...\")\n",
    "vcf_samples = [sample for sample in vcf.header.samples]\n",
    "vcf_egan_ids_with_rna = set(egan_ids_with_rna).intersection(set(vcf_samples))\n",
    "vcf.subset_samples(vcf_egan_ids_with_rna)\n",
    "vcf_samples_with_rna = [sample for sample in vcf.header.samples]\n",
    "print(f\"Number of samples in VCF with RNA ID and passing QC: {len(vcf_samples_with_rna)}\")\n",
    "# subset egan ID and RNA ID links to samples with SV calls and passing QC \n",
    "wgs_rna_paired_smpls_with_sv_calls_df = wgs_rna_paired_smpls_df[wgs_rna_paired_smpls_df.egan_id.isin(vcf_samples_with_rna)]\n",
    "# write EGAN-RNA ID pairs to file\n",
    "egan_rna_smpls_dir = root_dir_path.joinpath(\"egan_rna_smpls\")\n",
    "Path(egan_rna_smpls_dir).mkdir(parents=True, exist_ok=True)\n",
    "wgs_rna_paired_smpls_with_sv_calls_df.to_csv(egan_rna_smpls_dir.joinpath(\"egan_rna_ids_paired_pass_qc.tsv\"), sep=\"\\t\", index=False)\n",
    "rna_id_pass_qc_sv_calls = wgs_rna_paired_smpls_with_sv_calls_df.rna_id.unique().tolist()\n",
    "print(f\"Number of RNA IDs passing QC: {len(rna_id_pass_qc_sv_calls)}\")\n",
    "\n",
    "### write bed file for genes on chromosome passing QC \n",
    "gene_bed_dir = root_dir_path.joinpath(\"genes_bed\")\n",
    "gene_bed_dir.mkdir(parents=True, exist_ok=True)\n",
    "gene_bed_path = gene_bed_dir.joinpath(f\"{chrom}_genes.bed\")\n",
    "gene_id_pass_qc_on_chrom = []\n",
    "with open(gene_bed_path, \"w\") as f:\n",
    "    for gtf in pysam.TabixFile(gencode_path).fetch(chrom, parser = pysam.asGTF()):\n",
    "        if gtf.feature == \"gene\" and gtf.gene_id.split('.')[0] in gene_id_pass_qc_set:\n",
    "            # check for multiple entries with same name \n",
    "            gene_id_list, chrom_list, start_list, end_list, strand_list = [], [], [], [], []\n",
    "            gene_id_list.append(gtf.gene_id.split('.')[0])\n",
    "            chrom_list.append(gtf.contig)\n",
    "            start_list.append(gtf.start)\n",
    "            end_list.append((gtf.end))\n",
    "            strand_list.append((gtf.strand))\n",
    "            # check or write output\n",
    "            if len(gene_id_list) > 1 or len(chrom_list) > 1 or len(start_list) > 1 or len(end_list) > 1:\n",
    "                print(f\"{gene_id} has multiple entries in gencode file - excluded from output file.\")\n",
    "            elif len(chrom_list) == 0 or len(start_list) == 0 or len(end_list) == 0: \n",
    "                print(f\"{gene_id} has no entries in gencode file - excluded from output file.\")\n",
    "            else: \n",
    "                gene_id, gtf_chrom, start, end, strand = gene_id_list[0], chrom_list[0], start_list[0], end_list[0], strand_list[0]\n",
    "                gene_id_pass_qc_on_chrom.append(gene_id)\n",
    "                chrom_num = gtf_chrom.split(\"chr\")[1]\n",
    "                f.write(f\"{chrom_num}\\t{start}\\t{end}\\t{gene_id}\\t0\\t{strand}\\n\")\n",
    "print(f\"Number of genes pass QC on {chrom}: {len(gene_id_pass_qc_on_chrom)}\")\n",
    "\n",
    "# subset gene expression file to genes on chromosome \n",
    "ge_matrix_flat_chrom_df = ge_matrix_flat_df[ge_matrix_flat_df.gene_id.isin(gene_id_pass_qc_on_chrom)]\n",
    "# subset gene expression file to samples with SV calls \n",
    "ge_matrix_flat_chrom_egan_df = pd.merge(ge_matrix_flat_chrom_df, wgs_rna_paired_smpls_with_sv_calls_df, how=\"inner\", on=\"rna_id\")\n",
    "print(f\"Sample IDs in gene expression matrix with SV calls: {len(ge_matrix_flat_chrom_egan_df.egan_id.unique())}\")\n",
    "# write to file \n",
    "ge_matrix_flat_chrom_egan_dir = root_dir_path.joinpath(\"express_mtx\")\n",
    "Path(ge_matrix_flat_chrom_egan_dir).mkdir(parents=True, exist_ok=True)\n",
    "ge_matrix_flat_chrom_egan_df.to_csv(ge_matrix_flat_chrom_egan_dir.joinpath(f\"{chrom}_ge_mtx_flat.tsv\"), sep=\"\\t\", index=False)\n",
    "# add gene-sample pair to expression dataframe\n",
    "ge_matrix_flat_chrom_egan_df[\"gene_smpl_pair\"] = ge_matrix_flat_chrom_egan_df.gene_id + \",\" + ge_matrix_flat_chrom_egan_df.rna_id\n",
    "\n",
    "### write SVs on chromosome to bed file \n",
    "vrnts_bed_dir = root_dir_path.joinpath(\"vrnts_bed\")\n",
    "vrnts_bed_dir.mkdir(parents=True, exist_ok=True)\n",
    "vrnts_bed_path = vrnts_bed_dir.joinpath(f\"{chrom}_vrnts.bed\")\n",
    "with open(sv_info_path, \"r\") as f_in, open(vrnts_bed_path, \"w\") as f_out:\n",
    "    for line in f_in:\n",
    "        if line.startswith(\"plinkID\"): \n",
    "            continue\n",
    "        else: \n",
    "            vrnt_id, sv_chrom, pos, end = line.split(\"\\t\")[0], line.split(\"\\t\")[2], line.split(\"\\t\")[3], line.split(\"\\t\")[4]\n",
    "            if sv_chrom == chrom:\n",
    "                chrom_num = sv_chrom.split(\"chr\")[1]\n",
    "                f_out.write(f\"{chrom_num}\\t{pos}\\t{end}\\t{vrnt_id}\\n\")\n",
    "\n",
    "### SV information \n",
    "sv_info_df = pd.read_csv(sv_info_path, sep=\"\\t\", dtype={\"plinkID\":str})\n",
    "sv_types_list = sv_info_df.SVTYPE.unique().tolist()\n",
    "print(f\"SV types in SV info file: {sv_types_list}\")\n",
    "sv_info_id_af_df = sv_info_df[[\"plinkID\", \"AF\", \"SVTYPE\"]].rename(columns={\"plinkID\":\"vrnt_id\"})\n",
    "\n",
    "### generate output directories \n",
    "# variant-window intersection directory \n",
    "intersect_bed_dir=root_dir_path.joinpath(\"intersect_bed\")\n",
    "Path(intersect_bed_dir).mkdir(parents=True, exist_ok=True)\n",
    "# genotypes of variants in windows directory \n",
    "intersect_vrnt_gts_dir = root_dir_path.joinpath(\"intersect_vrnt_gts\")\n",
    "Path(intersect_vrnt_gts_dir).mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "# load variant and gene bed files \n",
    "vrnts_bed = BedTool(vrnts_bed_path)\n",
    "genes_bed = BedTool(gene_bed_path)\n",
    "# column names for variant window intersection bed file \n",
    "intersect_bed_columns={0:\"chrom_gene\", 1:\"start_gene\", 2:\"end_gene\", 3: \"gene_id\", 4: \"score\", 5: \"strand\", \n",
    "                       6:\"chrom_vrnt\", 7:\"start_vrnt\", 8:\"end_vrnt\", 9:\"vrnt_id\"}\n",
    "\n",
    "# load most severe consequence \n",
    "vep_msc_df = pd.read_csv(vep_msc_path, sep=\"\\t\", dtype={\"vrnt_id\": str})\n",
    "msc_list = vep_msc_df.Consequence.unique().tolist()\n",
    "vep_msc_cnsqn_df = vep_msc_df.rename(columns={\"Uploaded_variation\": \"vrnt_id\"})[[\"vrnt_id\", \"Consequence\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1b7922bf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>21</th>\n",
       "      <th>46151613</th>\n",
       "      <th>6374552</th>\n",
       "      <th>241388</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>7920738</td>\n",
       "      <td>241401</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>7925482</td>\n",
       "      <td>241407</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>7939915</td>\n",
       "      <td>241411</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>9737053</td>\n",
       "      <td>241454</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>9759937</td>\n",
       "      <td>241457</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1465</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>46618369</td>\n",
       "      <td>DEL_chr21_46616504_46618369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1466</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>46632239</td>\n",
       "      <td>DEL_chr21_46629641_46632239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1467</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>46668277</td>\n",
       "      <td>DEL_chr21_46637066_46668277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1468</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>9838378</td>\n",
       "      <td>DEL_chr21_9836915_9838378</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1469</th>\n",
       "      <td>21</td>\n",
       "      <td>46151613</td>\n",
       "      <td>10050232</td>\n",
       "      <td>DEL_chr21_9994239_10050232</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1470 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      21  46151613   6374552                       241388\n",
       "0     21  46151613   7920738                       241401\n",
       "1     21  46151613   7925482                       241407\n",
       "2     21  46151613   7939915                       241411\n",
       "3     21  46151613   9737053                       241454\n",
       "4     21  46151613   9759937                       241457\n",
       "...   ..       ...       ...                          ...\n",
       "1465  21  46151613  46618369  DEL_chr21_46616504_46618369\n",
       "1466  21  46151613  46632239  DEL_chr21_46629641_46632239\n",
       "1467  21  46151613  46668277  DEL_chr21_46637066_46668277\n",
       "1468  21  46151613   9838378    DEL_chr21_9836915_9838378\n",
       "1469  21  46151613  10050232   DEL_chr21_9994239_10050232\n",
       "\n",
       "[1470 rows x 4 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.read_csv(vrnts_bed_path, sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6e979cd1",
   "metadata": {},
   "outputs": [
    {
     "ename": "BEDToolsError",
     "evalue": "\nCommand was:\n\n\tbedtools window -sw -r 0 -l 0 -b /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/vrnts_bed/chr21_vrnts.bed -a /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/genes_bed/chr21_genes.bed\n\nError message was:\nError: malformed BED entry at line 1. Start was greater than end. Exiting.\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mBEDToolsError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_97886/4049522318.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m                                                            \u001b[0ml\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0ml\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m                                                            \u001b[0mr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m                                                            \u001b[0msw\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m                                                           ))) \n",
      "\u001b[0;32m~/.local/lib/python3.7/site-packages/pybedtools/bedtool.py\u001b[0m in \u001b[0;36mdecorated\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    915\u001b[0m             \u001b[0;31m# this calls the actual method in the first place; *result* is\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    916\u001b[0m             \u001b[0;31m# whatever you get back\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 917\u001b[0;31m             \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    918\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    919\u001b[0m             \u001b[0;31m# add appropriate tags\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.7/site-packages/pybedtools/bedtool.py\u001b[0m in \u001b[0;36mwrapped\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    399\u001b[0m                 \u001b[0mstdin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstdin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    400\u001b[0m                 \u001b[0mcheck_stderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcheck_stderr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 401\u001b[0;31m                 \u001b[0mdecode_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdecode_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    402\u001b[0m             )\n\u001b[1;32m    403\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.7/site-packages/pybedtools/helpers.py\u001b[0m in \u001b[0;36mcall_bedtools\u001b[0;34m(cmds, tmpfn, stdin, check_stderr, decode_output, encode_input)\u001b[0m\n\u001b[1;32m    458\u001b[0m                 \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwrite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    459\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 460\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mBEDToolsError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msubprocess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlist2cmdline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcmds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    461\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    462\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mOSError\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mIOError\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0merr\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mBEDToolsError\u001b[0m: \nCommand was:\n\n\tbedtools window -sw -r 0 -l 0 -b /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/vrnts_bed/chr21_vrnts.bed -a /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/genes_bed/chr21_genes.bed\n\nError message was:\nError: malformed BED entry at line 1. Start was greater than end. Exiting.\n"
     ]
    }
   ],
   "source": [
    "vrnts_window_intersect_str = StringIO(str(genes_bed.window(vrnts_bed, \n",
    "                                                           l=l, \n",
    "                                                           r=r,\n",
    "                                                           sw=True\n",
    "                                                          ))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8cdeeba9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Counting variants upstream, 0bp ...\n"
     ]
    },
    {
     "ename": "BEDToolsError",
     "evalue": "\nCommand was:\n\n\tbedtools window -sw -r 0 -l 0 -b /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/vrnts_bed/chr21_vrnts.bed -a /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/genes_bed/chr21_genes.bed\n\nError message was:\nError: malformed BED entry at line 1. Start was greater than end. Exiting.\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mBEDToolsError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_97886/1987748182.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     13\u001b[0m                                                                    \u001b[0ml\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0ml\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m                                                                    \u001b[0mr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m                                                                    \u001b[0msw\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     16\u001b[0m                                                                   )))        \n\u001b[1;32m     17\u001b[0m         \u001b[0mintersect_bed_df\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvrnts_window_intersect_str\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"\\t\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrename\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mintersect_bed_columns\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.7/site-packages/pybedtools/bedtool.py\u001b[0m in \u001b[0;36mdecorated\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    915\u001b[0m             \u001b[0;31m# this calls the actual method in the first place; *result* is\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    916\u001b[0m             \u001b[0;31m# whatever you get back\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 917\u001b[0;31m             \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    918\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    919\u001b[0m             \u001b[0;31m# add appropriate tags\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.7/site-packages/pybedtools/bedtool.py\u001b[0m in \u001b[0;36mwrapped\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    399\u001b[0m                 \u001b[0mstdin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstdin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    400\u001b[0m                 \u001b[0mcheck_stderr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcheck_stderr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 401\u001b[0;31m                 \u001b[0mdecode_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdecode_output\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    402\u001b[0m             )\n\u001b[1;32m    403\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.7/site-packages/pybedtools/helpers.py\u001b[0m in \u001b[0;36mcall_bedtools\u001b[0;34m(cmds, tmpfn, stdin, check_stderr, decode_output, encode_input)\u001b[0m\n\u001b[1;32m    458\u001b[0m                 \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwrite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstderr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    459\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 460\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mBEDToolsError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msubprocess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlist2cmdline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcmds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    461\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    462\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mOSError\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mIOError\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0merr\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mBEDToolsError\u001b[0m: \nCommand was:\n\n\tbedtools window -sw -r 0 -l 0 -b /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/vrnts_bed/chr21_vrnts.bed -a /lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/4_vrnt_enrich/sv_count_carriers/test/windows/200kb_window/genes_bed/chr21_genes.bed\n\nError message was:\nError: malformed BED entry at line 1. Start was greater than end. Exiting.\n"
     ]
    }
   ],
   "source": [
    "entry = 0 \n",
    "carrier_count = {}\n",
    "for direction in [\"upstream\", \"downstream\"]:\n",
    "    window_size = window_start\n",
    "    seen_sv_gene_pairs = set()\n",
    "    while window_size <= window_max:\n",
    "        print(f\"Counting variants {direction}, {window_size}bp ...\")\n",
    "        if direction == \"upstream\": \n",
    "            l, r = window_size, 0 \n",
    "        else: \n",
    "            l, r = 0, window_size\n",
    "        vrnts_window_intersect_str = StringIO(str(genes_bed.window(vrnts_bed, \n",
    "                                                                   l=l, \n",
    "                                                                   r=r,\n",
    "                                                                   sw=True\n",
    "                                                                  )))        \n",
    "        intersect_bed_df = pd.read_csv(vrnts_window_intersect_str, sep=\"\\t\", header=None).rename(columns=intersect_bed_columns)\n",
    "        # variant gene ID column \n",
    "        intersect_bed_df[\"vrnt_gene_id\"] = intersect_bed_df[\"vrnt_id\"].astype(str) + \",\" + intersect_bed_df[\"gene_id\"].astype(str)\n",
    "        # subset to variant gene pairs that have not been seen in previous windows    \n",
    "        sv_gene_pairs_in_window = set(intersect_bed_df.vrnt_gene_id.unique())\n",
    "        new_sv_gene_pairs = sv_gene_pairs_in_window - seen_sv_gene_pairs\n",
    "        if new_sv_gene_pairs.union(seen_sv_gene_pairs) != sv_gene_pairs_in_window: \n",
    "            raise ValueError(\"Missing variants from windows.\")\n",
    "        intersect_bed_new_sv_gene_df = intersect_bed_df[intersect_bed_df[\"vrnt_gene_id\"].isin(new_sv_gene_pairs)] \n",
    "        # generate set of variant IDs inside gene windows\n",
    "        intersect_vrnt_ids_no_dupl_df = intersect_bed_new_sv_gene_df[[\"chrom_vrnt\", \"vrnt_id\", \"start_vrnt\", \"end_vrnt\"]].drop_duplicates()\n",
    "        intersect_vrnt_ids_list = intersect_vrnt_ids_no_dupl_df.vrnt_id.unique()\n",
    "        # genotypes for variant IDs in windows\n",
    "        vrnt_gt_egan_dict = {}\n",
    "        count = 0\n",
    "        for vrnt_id in intersect_vrnt_ids_list:\n",
    "            # get chromosome, start and end of SV \n",
    "            chrom_vrnt, pos, end, = [intersect_vrnt_ids_no_dupl_df[intersect_vrnt_ids_no_dupl_df.vrnt_id == vrnt_id][col].item() for col in [\"chrom_vrnt\", \"start_vrnt\", \"end_vrnt\"]]\n",
    "            records = vcf.fetch(str(chrom_vrnt), pos, end)\n",
    "            found_vrnt_id = False\n",
    "            for rec in records: \n",
    "                vcf_vrnt_id = str(rec.id)\n",
    "                if vrnt_id == vcf_vrnt_id:\n",
    "                    found_vrnt_id = True \n",
    "                    gts = [s[\"GT\"] for s in rec.samples.values()]\n",
    "                    for i, gt in enumerate(gts): \n",
    "                        vrnt_gt_egan_dict[count] = [vrnt_id, vcf.header.samples[i], gt]\n",
    "                        count += 1 \n",
    "            if not found_vrnt_id: \n",
    "                raise ValueError(f\"Did not find {vrnt_id} in {vcf_path}\")\n",
    "        vrnt_gt_egan_nogene_df = pd.DataFrame.from_dict(vrnt_gt_egan_dict, orient=\"index\", columns=[\"vrnt_id\", \"egan_id\", \"genotype\"])\n",
    "        vrnt_gt_egan_nogene_df = vrnt_gt_egan_nogene_df.astype({\"genotype\":str})\n",
    "\n",
    "        # merge \n",
    "        dfs_to_merge = [vrnt_gt_egan_nogene_df, \n",
    "                        intersect_bed_new_sv_gene_df[[\"vrnt_id\", \"gene_id\",\"vrnt_gene_id\"]].drop_duplicates(),\n",
    "                        sv_info_id_af_df, \n",
    "                        vep_msc_cnsqn_df\n",
    "                       ]\n",
    "        df_merged = reduce(lambda  left,right: pd.merge(left,right,on=['vrnt_id'],\n",
    "                                                    how='inner'), dfs_to_merge)\n",
    "        # add expression \n",
    "        sv_intersect_express_info_df = pd.merge(df_merged, \n",
    "                                             ge_matrix_flat_chrom_egan_df,                  \n",
    "                                             how=\"inner\",\n",
    "                                             on=[\"egan_id\", \"gene_id\"])\n",
    "        # write to file \n",
    "        intersect_vrnt_gts_path = intersect_vrnt_gts_dir.joinpath(f\"{chrom}_intersect_vrnt_gts_{direction}_{window_size}.tsv\")\n",
    "        sv_intersect_express_info_df.to_csv(intersect_vrnt_gts_path, sep=\"\\t\", index=False)\n",
    "\n",
    "        ### count carriers in test and controls across different Z-score cutoffs \n",
    "        for z_cutoff in z_cutoff_list: \n",
    "            carrier_count[entry] = [chrom, direction, window_size, z_cutoff]\n",
    "            ### get control and test gene-sample pairs based on TPM cutoff \n",
    "\n",
    "            # get gene-sample pairs passing z-score and TPM cutoff \n",
    "            misexp_df = ge_matrix_flat_chrom_egan_df[(ge_matrix_flat_chrom_egan_df[\"z-score\"] > z_cutoff) &\n",
    "                                                     (ge_matrix_flat_chrom_egan_df[\"TPM\"] > tpm_cutoff)\n",
    "                                                    ]\n",
    "            test_gene_ids = misexp_df.gene_id.unique()\n",
    "            test_gene_smpl_pairs = misexp_df.gene_smpl_pair.unique()\n",
    "            # get gene-sample pairs in control group - same gene set, samples not misexpressing\n",
    "            cntrl_gene_smpl_pairs = ge_matrix_flat_chrom_egan_df[~ge_matrix_flat_chrom_egan_df.gene_smpl_pair.isin(test_gene_smpl_pairs) &\n",
    "                                                                 ge_matrix_flat_chrom_egan_df.gene_id.isin(test_gene_ids)\n",
    "                                                                ].gene_smpl_pair.unique()\n",
    "            # check overlap between sets is empty \n",
    "            if len(set(test_gene_smpl_pairs).intersection(set(cntrl_gene_smpl_pairs))) != 0:\n",
    "                raise ValueError(\"Overlap between control and test gene-sample pair sets\")\n",
    "            # check number of test and control gene-pairs matches expected \n",
    "            total_test_control_pairs = len(test_gene_smpl_pairs) + len(cntrl_gene_smpl_pairs)\n",
    "            test_genes_by_smpls = len(test_gene_ids) * len(rna_id_pass_qc_sv_calls)\n",
    "            if total_test_control_pairs !=  test_genes_by_smpls: \n",
    "                raise ValueError(\"Number of test and control gene pairs does not match test gene IDs by samples\")\n",
    "            carrier_count[entry] += [len(test_gene_ids), len(rna_id_pass_qc_sv_calls), len(test_gene_smpl_pairs), len(cntrl_gene_smpl_pairs)]\n",
    "\n",
    "            ### count carriers in control and test groups \n",
    "            # create column with gene-sample pairs  \n",
    "            sv_intersect_express_info_df[\"gene_smpls_id\"] = sv_intersect_express_info_df.gene_id.astype(str) + \",\" + sv_intersect_express_info_df.rna_id.astype(str)\n",
    "            # get misexpression carriers \n",
    "            misexpress_carriers_df = sv_intersect_express_info_df[(sv_intersect_express_info_df.gene_smpls_id.isin(test_gene_smpl_pairs)) & \n",
    "                                                                  (sv_intersect_express_info_df.AF >= af_lower) & \n",
    "                                                                  (sv_intersect_express_info_df.AF < af_upper) &\n",
    "                                                                  (sv_intersect_express_info_df.genotype.isin(['(0, 1)', '(1, 1)']))\n",
    "                                                                 ].copy()\n",
    "            # get non-misexpression carriers\n",
    "            cntrl_carriers_df = sv_intersect_express_info_df[(sv_intersect_express_info_df.gene_smpls_id.isin(cntrl_gene_smpl_pairs)) &\n",
    "                                                                      (sv_intersect_express_info_df.AF >= af_lower) & \n",
    "                                                                      (sv_intersect_express_info_df.AF < af_upper) &\n",
    "                                                                      (sv_intersect_express_info_df.genotype.isin(['(0, 1)', '(1, 1)']))\n",
    "                                                                     ].copy()\n",
    "            # add number of carriers \n",
    "            gene_smpl_misexpress_carriers = misexpress_carriers_df.gene_smpls_id.unique()\n",
    "            gene_smpl_cntrl_carriers = cntrl_carriers_df.gene_smpls_id.unique()\n",
    "            carrier_count[entry] += [len(gene_smpl_misexpress_carriers), len(gene_smpl_cntrl_carriers)]\n",
    "            # check overlap between sets is empty \n",
    "            if len(set(gene_smpl_misexpress_carriers).intersection(set(gene_smpl_cntrl_carriers))) != 0:\n",
    "                raise ValueError(\"Overlap between control and test gene-sample pair carrier sets\")\n",
    "\n",
    "            ### count carriers for SV types and most severe consequence \n",
    "            # question remains whether to assign priority to different SV classes\n",
    "            for sv_type in sv_types_list:\n",
    "                misexpress_carriers_sv_type_df = misexpress_carriers_df[misexpress_carriers_df.SVTYPE == sv_type]\n",
    "                cntrl_carriers_sv_type_df = cntrl_carriers_df[cntrl_carriers_df.SVTYPE == sv_type]\n",
    "                sv_type_misexp_carriers = len(misexpress_carriers_sv_type_df.gene_smpls_id.unique())\n",
    "                sv_type_cntrl_carriers = len(cntrl_carriers_sv_type_df.gene_smpls_id.unique())\n",
    "                carrier_count[entry] += [sv_type_misexp_carriers, sv_type_cntrl_carriers]\n",
    "                for msc in msc_list: \n",
    "                    msc_type_misexp_carriers = len(misexpress_carriers_sv_type_df[misexpress_carriers_sv_type_df.Consequence == msc].gene_smpls_id.unique())\n",
    "                    msc_type_cntrl_carriers = len(cntrl_carriers_sv_type_df[cntrl_carriers_sv_type_df.Consequence == msc].gene_smpls_id.unique())\n",
    "                    carrier_count[entry] += [msc_type_misexp_carriers, msc_type_cntrl_carriers]\n",
    "            entry += 1\n",
    "        # add variant gene pairs to seen set \n",
    "        seen_sv_gene_pairs = seen_sv_gene_pairs.union(sv_gene_pairs_in_window)\n",
    "        window_size += window_step_size\n",
    "\n",
    "# write results \n",
    "vrnt_type_cols = []\n",
    "for sv_type in sv_types_list: \n",
    "    vrnt_type_cols += [f\"{sv_type}_misexp\", f\"{sv_type}_contrl\"]\n",
    "    for msc in msc_list:\n",
    "        vrnt_type_cols += [f\"{sv_type}_{msc}_misexp\", f\"{sv_type}_{msc}_contrl\"]\n",
    "\n",
    "columns = [\"chrom\", \"direction\", \"window_size\", \"z_cutoff\", \"misexp_genes\", \"smpls_pass_qc\", \"total_misexp\", \"total_control\", \"all_sv_misexp\", \"all_sv_contrl\"] + vrnt_type_cols\n",
    "carrier_count_df = pd.DataFrame.from_dict(carrier_count, orient=\"index\", columns=columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62975f54",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
